Basic spatial data

March 2018, June 2019. Mauro Alberti.

Last run: 2019-06-16

Developement code:



In [1]:

    
%load_ext autoreload
%autoreload 1

1. Introduction

gsf is a library for the processing of geometric and geographic data, with a focus on structural geology data. It is composed by a Python 3 module, pygsf, and an experimental, still in progress Haskell module, named hsgsf. This notebook will present the Python version.

Since we will plot geometric data into stereonets, prior to any other operation, we import mplstereonet and run the IPython command %matplotlib inline, that allows to incorporate Matplotlib plots in the Jupyter notebook.



In [2]:

    
%matplotlib inline

We import all classes/methods from the geometry sub-module:



In [3]:

    
from pygsf.spatial.vectorial.geometries import *

2. Basic spatial data types: points and planes

The reference axis orientations used in pygsf are the x axis parallel to East, y parallel to North and z vertical, upward-directed.

2.1 Cartesian points

A point is created by providing three Cartesian coordinates:



In [4]:

    
p1 = Point(1.0, 2.4, 0.2)  # definition of a Point instance

When considering two points, we can calculate their 3D distance as well as their horizontal, 2D distance:



In [5]:

    
p2 = Point(0.9, 4.2, 10.5)



In [6]:

    
p1.dist3DWith(p2)  # 3D distance between two points



In [7]:

    
p1.dist2DWith(p2)  # horizontal (2D) distance between two points

Among other methods, we can:

translate the point position by providing three offset cartesian values (x, y and z) or directly via a vector;
check if two points are within a given range of each other;
convert a point to a vector.

2.2 Cartesian planes

A Cartesian plane can be defined in a few different ways, with the simplest one by providing three points within the plane:



In [8]:

    
pl1 = CPlane.fromPoints(Point(0, 0, 0, epsg_cd = 2000), Point(1, 0, 0, epsg_cd = 2000), Point(0, 1, 0, epsg_cd = 2000))



In [9]:

    
print(pl1)









    



CPlane(0.0000, 0.0000, 1.0000, 0.0000, 2000)

The four coefficient returned (a, b, c and d) define the Cartesian plane as in the equation:

ax + by + cz = d

For the given example, the equation is satisfied for all x and y values provided z is zero. We are therefore considering a horizontal plane passing through the frame origin.

The versor normal to a Cartesian plane is obtained by the method:



In [10]:

    
normal_versor = pl1.normVersor()  # versor (unit vector) normal to the provided Cartesian plane



In [11]:

    
print(normal_versor)









    



Vect(0.0000, 0.0000, 1.0000, EPSG: 2000)

In this example the normal versor is vertical.

We can calculate the intersection, expressed as a versor, between two Cartesian planes:



In [12]:

    
pl1, pl2 = CPlane(1, 0, 0, 0, epsg_cd = 2000), CPlane(0, 0, 1, 0, epsg_cd = 2000)
inters_v = pl1.intersVersor(pl2)  # intersection versor between two Cartesian planes 
print(inters_v)









    



Vect(0.0000, -1.0000, 0.0000, EPSG: 2000)